Quality metrics per cell

colnames(colData(sce))

##  [1] "tissue"                                        
##  [2] "group"                                         
##  [3] "mRNA"                                          
##  [4] "well"                                          
##  [5] "sex"                                           
##  [6] "age"                                           
##  [7] "diameter"                                      
##  [8] "cellId"                                        
##  [9] "level1class"                                   
## [10] "level2class"                                   
## [11] "is_cell_control"                               
## [12] "total_features_by_counts"                      
## [13] "log10_total_features_by_counts"                
## [14] "total_counts"                                  
## [15] "log10_total_counts"                            
## [16] "pct_counts_in_top_50_features"                 
## [17] "pct_counts_in_top_100_features"                
## [18] "pct_counts_in_top_200_features"                
## [19] "pct_counts_in_top_500_features"                
## [20] "total_features_by_counts_endogenous"           
## [21] "log10_total_features_by_counts_endogenous"     
## [22] "total_counts_endogenous"                       
## [23] "log10_total_counts_endogenous"                 
## [24] "pct_counts_endogenous"                         
## [25] "pct_counts_in_top_50_features_endogenous"      
## [26] "pct_counts_in_top_100_features_endogenous"     
## [27] "pct_counts_in_top_200_features_endogenous"     
## [28] "pct_counts_in_top_500_features_endogenous"     
## [29] "total_features_by_counts_feature_control"      
## [30] "log10_total_features_by_counts_feature_control"
## [31] "total_counts_feature_control"                  
## [32] "log10_total_counts_feature_control"            
## [33] "pct_counts_feature_control"                    
## [34] "pct_counts_in_top_50_features_feature_control" 
## [35] "pct_counts_in_top_100_features_feature_control"
## [36] "pct_counts_in_top_200_features_feature_control"
## [37] "pct_counts_in_top_500_features_feature_control"
## [38] "total_features_by_counts_Spike"                
## [39] "log10_total_features_by_counts_Spike"          
## [40] "total_counts_Spike"                            
## [41] "log10_total_counts_Spike"                      
## [42] "pct_counts_Spike"                              
## [43] "pct_counts_in_top_50_features_Spike"           
## [44] "pct_counts_in_top_100_features_Spike"          
## [45] "pct_counts_in_top_200_features_Spike"          
## [46] "pct_counts_in_top_500_features_Spike"          
## [47] "total_features_by_counts_Mt"                   
## [48] "log10_total_features_by_counts_Mt"             
## [49] "total_counts_Mt"                               
## [50] "log10_total_counts_Mt"                         
## [51] "pct_counts_Mt"                                 
## [52] "pct_counts_in_top_50_features_Mt"              
## [53] "pct_counts_in_top_100_features_Mt"             
## [54] "pct_counts_in_top_200_features_Mt"             
## [55] "pct_counts_in_top_500_features_Mt"

colData(sce)$total_counts[1:5] #Total count for each cell

## [1] 29060 29304 38929 40557 27625

colData(sce)$total_features_by_counts[1:5] #Number of genes with non zero counts per cell

## [1] 4898 4750 6090 5887 4753

if(length(is.spike==TRUE)>0){colData(sce)$pct_counts_Spike[1:5]} #Percentage of spike-in counts per cell

## [1] 23.07639 21.95946 16.27322 17.33856 21.46968

Quality metrics per gene

colnames(rowData(sce))

## [1] "is_feature_control"       "is_feature_control_Spike"
## [3] "is_feature_control_Mt"    "mean_counts"             
## [5] "log10_mean_counts"        "n_cells_by_counts"       
## [7] "pct_dropout_by_counts"    "total_counts"            
## [9] "log10_total_counts"

rowData(sce)$mean_counts[1:5] #Average count per gene

## [1] 0.27886855 0.42362729 1.04292845 0.04459235 0.37936772

rowData(sce)$n_cells_by_counts[1:5] #Number of cells with non zero counts per gene

## [1]  472  573 1167   77  669

rowData(sce)$pct_dropout_by_counts[1:5] #Percentage of dropouts per gene

## [1] 84.29285 80.93178 61.16473 97.43760 77.73710

Plot quality metrics

One can plot histograms based on quality metrics to determine possible thresholds for cell/gene filterings. We can also visualize some of the quality metrics on the PCA projected data.

assay(sce, "logcounts_raw") <- log2(counts(sce)+1) 

hist(sce$total_counts, breaks=50,  xlab="Library size", ylab="Number of cells")

plotPCA(sce, run_args=c(exprs_values="logcounts_raw"), size_by="total_counts")

hist(sce$total_features_by_counts, xlab="Number of expressed genes", breaks=50,  ylab="Number of cells")

plotPCA(sce, run_args=c(exprs_values="logcounts_raw"), size_by="total_features_by_counts")

if(length(is.spike==TRUE)>0)
{hist(sce$pct_counts_Spike, xlab="ERCC proportion (%)", ylab="Number of cells", breaks=50)
plotPCA(sce, run_args=c(exprs_values="logcounts_raw"), size_by="pct_counts_Spike")}

Automatic filtering of outlier cells based on quality metrics

Alternatively, isOutlier function provides an automatic detection of low quality cells based on quality metrics. In this example, we are filtering cells based on library size, number of expressed genes per cell and total count over all spike-in transcripts in each cell. Excluded cells are those with the total number of expressed genes and the total sum of counts more than 3 median absolute deviations below the median across the genes. Note that we dont filter many cells in this step, it is likely that dataset was already quality controlled by the original study.

libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", log=TRUE)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", log=TRUE)
if(length(is.spike==TRUE)>0){spike.drop <- isOutlier(sce$pct_counts_Spike, nmads=3, type="higher")
ind <- libsize.drop | feature.drop | spike.drop}else{ind <- libsize.drop | feature.drop}

sce$removed=ind
plotPCA(sce, colour_by="removed", run_args=c(exprs_values="logcounts_raw"), shape_by="removed")

if(length(is.spike==TRUE)>0)
{sce <- sce[,!(libsize.drop | feature.drop | spike.drop)]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop), BySpike=sum(spike.drop), Remaining=ncol(sce))
}else{sce <- sce[,!(libsize.drop | feature.drop)]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop), Remaining=ncol(sce))
}

##   ByLibSize ByFeature BySpike Remaining
## 1         8         3       8      2989

Explore highly expressed genes

We can also explore highly and overexpressed genes. The overexpressed genes, are those genes which expression is higher than others by several magnitudes. Overexpressed genes can bias further procedures such as clustering, therefore, they should be removed before performing downstream analysis. Note that highly expressed genes mostly include spike-ins and mitochondrial genes. One of the most commonly overexpressed endogenous gene is “Malat1”.

plotHighestExprs(sce, exprs_values = "counts")

Filter overexpressed features

In this step we will filter spike-in and mitochondrial genes as well as “Malat1”.

if(length(is.spike==TRUE)>0){sce=sce[-which(grepl("^ERCC-", rownames(sce))),]}
sce=sce[-which(grepl("^mt-", rownames(sce))),]
sce <- sce[-which(rownames(sce)=="Malat1"),]
dim(sce)

## [1] 19804  2989

Filter lowly expressed genes based on the average count per gene

Now, we should filter out lowly expressed genes as they do not provide any insights into the underlying biology. In the filtering we removed lowly expressed genes that are genes with average expression count (adjusted by library size) equal to 0. Note that we dont filter many genes in this step, it is likely that dataset was already quality controlled by the original study.

rowData(sce)$ave.counts <- calcAverage(sce, exprs_values = "counts", use_size_factors=FALSE)
to.keep <- rowData(sce)$ave.counts > 0
sce <- sce[to.keep,]
summary(to.keep)

##    Mode   FALSE    TRUE 
## logical       2   19802

dim(sce)

## [1] 19802  2989

Perform standard normalization by computing counts per million (CPM)

After quality control and filtering, one should normalize the dataset to remove potential technical bias. One strategy is to use CPM (Count Per Million) normalization which is a correction to remove the noise related to sequencing depth. CPM divides each count by its total sum (across all the genes) and multiplies by one million. In this way, each cell has the same total sum of the counts. We can plot PCA on the log-transformed CPM normalized counts and compare with previous PCA on the log raw counts.

cpm(sce) <- calculateCPM(sce, use_size_factors=FALSE)
cpm(sce)[1:4,1:4]

##                V3       V4      V5        V6
## Tspan12    0.0000  0.00000   0.000  94.75081
## Tshz1    144.5226 47.89501   0.000  63.16720
## Fnbp1l   144.5226 47.89501 197.336 126.33441
## Adamts15   0.0000  0.00000   0.000   0.00000

assay(sce, "logcounts") <- log2(cpm(sce)+1) 
plotPCA(sce,colour_by="level1class",run_args=c(exprs_values="logcounts"))

set.seed(100)
plotTSNE(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))

plotUMAP(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))

Perform normalization using decnvolution method

Another strategy, more suitable for scRNAseq data, is deconvolution method implemented in scran package. The deconvolution method normalizes data by cells-pooled size factors that account for dropout biases. More details about this normalization technique can be found in study https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0947-7. Similarly, we will plot PCA on the log scran normalized counts.

set.seed(100)
clusters <- quickCluster(sce, min.size=200, min.mean=0.1, method="igraph")
sce <- computeSumFactors(sce, cluster=clusters, min.mean=0.1)
sce <- scater::normalize(sce, return_log = FALSE)
normcounts(sce)[1:4,1:4]

##                V3        V4       V5       V6
## Tspan12  0.000000 0.0000000 0.000000 1.278561
## Tshz1    1.894734 0.6425685 0.000000 0.852374
## Fnbp1l   1.894734 0.6425685 2.582446 1.704748
## Adamts15 0.000000 0.0000000 0.000000 0.000000

assay(sce, "logcounts") <- log2(normcounts(sce)+1) 
plotPCA(sce,colour_by="level1class",run_args=c(exprs_values="logcounts"))

Feature selection - to retrieve most variable genes

Feature selection step is aimed at preserving biologically relevant information in the dataset and improving computational efficiency of downstream analyses. In most of the cases, it seeks for highly variable genes based on mean/variance relationship. M3Drop is one of the packages for automatic detection of most variable genes. It first calculates the mean and square coefficient of variation and fits the quadratic curve between the two variables. Then a chi-square test is used to find high variable genes which are significantly above the curve. We found CPM adjusted counts more appriopriate to use as input for this method rather than scran normalized values. Note that when using older versions of M3Drop package (<1.10.0) the output is an array of genes instead of gene table with p-values.

hvg_table <- BrenneckeGetVariableGenes(cpm(sce), fdr = 0.01, minBiolDisp = 0.5)

length(hvg_table$Gene)

## [1] 8183

sce_sub=sce[as.character(hvg_table$Gene),]

Dimension reduction - to project data into low-dimensional embedding

Dimension reduction can be used to visualize basic structure of the data or to reduce the number of features prior downstream analysis such as clustering. In contrast to the feature selection which extracts the most informative features, dimension reduction techniques projects the data into a new low-dimensional space that preserves the structure of the data. To reduce dataset dimensionality one can use previously mentioned PCA or more novel techniques such as tSNE or UMAP (all implemented in scater package). Note that projections should be applied to normalized and log-transformed count matrices. PCA is a deterministic approach while tSNE and UMAP are stochastic - for reproducibility of tSNE and UMAP projections one has to set the seed for generating random variables.

assay(sce_sub, "logcounts") <- log2(cpm(sce_sub)+1) 
plotPCA(sce_sub, colour_by="level1class", run_args=c(exprs_values = "logcounts"))  

set.seed(100)
plotTSNE(sce_sub, colour_by="level1class", run_args=c(exprs_values = "logcounts"))

plotUMAP(sce_sub, colour_by="level1class", run_args=c(exprs_values = "logcounts"))

Clusterize cells using SC3 package

To cluster cells we can use SC3 package. SC3 is a consensus clustering method for single-cell RNA-seq data. It first calculates distances between the cells using three metrics: Euclidean, Pearson, and Spearman. On each of the obtained distance matrix two transformations are applied: PCA and Laplacian graph. Then K-Means clustering technique is used to cluster the transformed distance matrices subject to the first d eigenvectors. In result, several individual clusterings are obtained which are further combined into a single consensus clustering using Cluster-based Similarity Partitioning Algorithm (CSPA). In this example we will use available cell annotation in the clustering inference of SC3. In the real life applications we usually do not know the exact number populations in the sample. Note that clustering cells with SC3 may take a substantial amount of time. We will visualize results on the PCA projected space.

k_input=length(unique(sce$level1class))
rowData(sce)$feature_symbol <- rownames(sce)
sce <- sc3(sce, ks = k_input, gene_filter = FALSE, biology = FALSE)

table(colData(sce)[paste0("sc3_",k_input, "_clusters")][[1]])

## 
##   1   2   3   4   5   6   7 
## 960   8 708  49 517 311 436

assay(sce, "logcounts") <- log2(cpm(sce)+1)
plotTSNE(sce, colour_by=paste0("sc3_",k_input, "_clusters"), run_args=c(exprs_values = "logcounts"))

adjustedRandIndex(colData(sce)[paste0("sc3_",k_input, "_clusters")][[1]],sce$level1class)

## [1] 0.7940186

Exercise

How the clustering performance would change if you would estimate the number of cell populations in SC3 using sc3_estimate_k function? Is the number of estimated clusters close to the annotated number of cell populations ?